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Abstract 



0^ ^Ve developed a Monte Carlo event generator for production of nucleon configurations in complex nuclei consistently 
including effects of Nucleon-Nucleon (NN) correlations. Our approach is based on the Metropolis search for configurations 



satisfying essential constraints imposed by short- and long-range NN correlations, guided by the findings of realistic 



C/3 



calculations of one- and two-body densities for medium-heavy nuclei. The produced event generator can be used for 
Monte Carlo (MC) studies of pA and AA collisions. We perform several tests of consistency of the code and comparison 
with previous models, in the case of high energy proton- nucleus scattering on an event-by-event basis, using nucleus 
configurations produced by our code and Glauber multiple scattering theory both for the uncorrelated and the correlated 
^ configurations; fluctuations of the average number of collisions are shown to be affected considerably by the introduction 

, , of NN correlations in the target nucleus. We also use the generator to estimate maximal possible gluon nuclear shadowing 

in a simple geometric model. 
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1. Introduction 

The structure of nuclei has been described for a long 
time by independent particle models in which the nu- 
cleus is treated as a collection of fermions freely moving 
about, this picture being suggested by the fact that nu- 
clear matter is a dilute system. Nevertheless the short 
range structure of nuclei cannot be accurately described by 
this simplified picture, and essential aspects of nuclei such 
as high momentum component of the nuclear wave func- 
tion are completely missed by the non-interacting model; 
configurations of nucleons at short separations are known 
as short range correlations (SRC) [J, Q|. NN correlations 
have been recently unambiguously observed^ in a series 
of dedicated experiments, in which high-momentum nu- 
cleons were knocked out by high energy probes and their 
correlated partner were observed with opposite momen- 
tum; the probability for a nucleon to belong to a SRC pair 
was measured to be about 20% in 12 C and 25% in in heavy 
nuclei; see Ref.Q for a review. A theoretical description 
of SRC can be made in a satisfactory way using nonrela- 
tivistic Hamiltonians containing two- and three-body real- 
istic potentials. It is currently possible to solve exactly the 
Schrodinger equation for light nuclei while, for A > 12, ap- 
proximate methods must be considered, which have been 
shown to produce a reasonable approximation to the basic 
ground state quantities such as one- and two-body densi- 
ties and momentum distributions. These studies demon- 
strated that NN correlations produce non-negligible ef- 
fects in a variety of processes including those which are 
not specially probing SRC. For example, significant ef- 



fects are found for high energy total nucleon-nucleus cross 
sections 0, and it has thus been shown that they must be 
taken into account in the theoretical calculation of such 
processes. In this Letter we are addressing the challenging 
problem of including NN correlations in the description 
of high energy proton-nucleus and nucleus- nucleus colli- 
sions. It was shown in Ref.Q that using correlated two- 
body densities for the analysis of fluctuations of the mean 
number of participating nucleons in proton-nucleus colli- 
sions significantly modifies the results, especially the tails 
of the distribution, and it was pointed out the importance 
of implementing NN correlations in the event generator 
for nucleon configurations in a way consistent with the 
single nucleon nuclear density and use it for description of 
the heavy ion collisions on an event-by-event basis. Sim- 
ulations of such processes are commonly performed using 
configurations of nucleons as an input, consisting of nu- 
cleons spatial locations and isospins, generated assuming 
the independent particle model; however, inclusion of NN 
correlations can provide a more realistic description of ini- 
tial states of nuclei and appears to be feasible with a good 
accuracy. 

The Letter is organized as follows; in section II, de- 
tails of NN correlations are discussed and the method 
adopted for generation of configurations is outlined; we 
describe the generation of nuclear configurations within 
a MC approach including central two-body correlations 
which are consistent with realistic calculations of one- and 
two-body density distributions. In section III we introduce 
the probability distribution functions Pn for the iV-th. or- 
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der, N = 1, A, interaction of a projectile nucleon with a 
nucleus A, and we present results for fluctuations. In sec- 
tion IV we use our event generator to calculate maximal 
possible shadowing for gluons based on simple geometric 
considerations. 

2. NN correlations 

The aim of the present work is to develop a MC proce- 
dure for generating spatial nucleon configurations to be 
used as an input for simulations of collisions involving 
heavy nuclei which treats NN correlations in a realistic 
way. The commonly used approach to this problem is tak- 
ing a given number of nucleons out of a Woods-Saxon den- 
sity distribution, which describes the probability density 
function of a nucleon to be located at a distance r from 
the center of the system. As a result, the A nucleons are 
positioned independently from each other, such that all 
two-particle (and higher) correlations are completely ig- 
nored; this is justified by the assumption that inclusive 
quantities are not very much dependent on these correla- 
tions. Alternatively one puts nucleons in consequently one 
after another imposing the condition that the distance be- 
tween nucleons should be larger than some minimal one. 
This procedure, however, results in a wrong single nucleon 
distribution of the nucleons and must be improved. 

We argue that when considering event-by-event observ- 
ables, NN correlations in the nuclear wave function are rel- 
evant. This extends observation of 6| that NN correlations 
significantly modify the variance of the distribution over 
the number of collisions. In the present work we analyze 
several different approximations for correlations. A full 
implementation of realistic short range correlations 0, Q 
in the wave function ^> of a complex nucleus can be made 
using an independent particle model wave function $ and 
a proper correlation operator F — nf<j f(fij) as follows: 



$>(n,...,r A ) = F<f>(r u ...,r A ). 



(1) 



/ being a set of state-dependent correlation functions 
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nonrelativistic Hamiltonian with realistic two- 



obtained variationally from a 
and three- 
body interactions; O^' is the spin-isospin operator ap- 
pearing in the realistic two-body potential, usually written 
in the form V = (>(r«) = £„ w<">(r y ) 6%] . 

In the present implementation of correlations, we did 
not consider the full realistic set of correlations shown in 
Eq-CJ, which would imply using the full \^f\ 2 as a prob- 
ability function for the Metropolis random search. This 
would be numerically very demanding, especially for heavy 
nuclei, as it would involve generating a full configuration 
of A nucleons. Hence as a first step we will implement 
only effective central correlations; the inclusion of spin- 
and isospin-dependent configurations, which allows to ac- 
count for most of the properties of one- and two-body den- 
sities described in 0, 0, [13] and references therein, will be 
discussed later on. 



A major problem with using a system of uncorrelated 
nucleons for the event generators is the occurrence of hot 
spots: regions in space where two or more nucleons are 
allowed to overlap and produce unrealistically high local 
densities. A simple excluded volume model, in which nu- 
cleons are prevented from being closer than a fixed distance 
of about the nucleon size d = 1 fm does not account for the 
realistic picture of nuclei, since this approach artificially 
moves density from the center of the nucleus towards the 
periphery, distorting the nuclear profile and altering the 
surface properties. The first constraint in our generator 
is thus to preserve the one body density of the nucleus, 
defined as follows: 



A 

A / JJ dr t \y(ri,...,r A ) 

J i=2 



(2) 



We suggest a method for generating the spatial configu- 
rations of nuclei with the following essential features: i) 
the nucleons arc distributed according to a given single- 
particle distribution; ii) each configuration embodies the 
NN pair correlations. The resulting two-body densities are 
by far more realistic than the commonly used independent 
particle model ones, or even those of the excluded volume 
approximation. The outlined results can be achieved by 
first generating a large number of random nucleon posi- 
tions in a cube with given density (0.17 nucleons/ fm 3 ). 
Next, a Metropolis random search is performed, with a 
probability function given by the square of the wave func- 
tion; at this stage, we impose the constraint of the Woods- 
Saxon probability distribution (Ref.ll!]), choosing those 
nucleons which are relevant to the desired one-body den- 
sity. The procedure comes to an end when exactly A nu- 
cleons satisfy the constraints at the same time. 

It was shown in Refs.0, Q that SRC have sizable effects 
on the single particle density (Eq.([2|)) as well as on those 
quantities which depend depend on the two-body density 
of the system, defined as 

r A 

p^(r l7 r 2 ) = A(A-l) / JJ dr t |*(n, r A )\ 2 , (3) 

^ «=3 

such as the two-body momentum distributions. In many 
theoretical studies, Eq.© is approximated as follows 



P (2) (ri,r 2 ) = 



A-l 
~A~ 



P (1) (n)p (1) (r 2 ) [1 - C(n,r 2 )] (4) 



where C(t*i, t* 2 ) vanishes if NN and statistical correlations 
are disregarded. We can compare the realistic results [71 [l2| 
for the radial two-body density, defined as the integral of 
Eq.Q over the pair center of mass, R: 

P {2) (r) = JdR (r ! = R + , r 2 = R - V) (5) 

and the results from our MC generator for the same quan- 
tity. If we assume that C(t*i, t* 2 ) depends only on the rela- 
tive coordinate r — \r\ — r 2 | (which holds true exactly for 
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Figure 1: The quantity of Eq.JgJl, calculated for le O. Dashed 
histogram: step function (excluded volume) correlation; solid his- 
togram: Gaussian correlation; solid curve: the realistic calculation 
of Ref. 0. 



infinite nuclear matter and which is a good approximation 
in finite nuclei, far enough from the edge of the nucleus), 
the pair distribution function C(r) can be obtained from 

(2) (2) 

the correlated p c (r) and uncorrelated p\j (r) two-body 
radial densities: 



C{r) = 1 



(6) 



The quantities defined in Eq.© and Eq.© were calcu- 
lated using the MC code using several approximations for 
the correlation function /(r) (see Eq.([T])). They are shown 
in Fig[IJ the results of Ref. [7( are also shown. Our proce- 
dure was to start with a central correlation function f(r) 
taken as a step function, i.e. vanishing wave function for 
r < 1 fm; we then used a Gaussian f(r) = 1 — exp(—f3r 2 ) 
in order to match better the realistic result. It is known 
that the realistic central correlations has a small overshoot- 
ing over unity, and it is in general different from our Gaus- 
sian correlation; however, in realistic calculations of finite 
nuclei 0, E3] the peak in the radial two-body density is 
mainly due to the existence of state-dependent, mainly 
tensor correlations, and three-body diagrams. For this 
reason we prefer not to introduce an unrealistic central 
correlation to reproduce such a peak and use a parameter 
(3 = 0.9 fm~ 2 in the Gaussian correlation which gives a 
healing distance similar to the one exhibited by the real- 
istic approach. This amounts to calculate the wave func- 
tion with a central correlation operator; the fully realistic 
implementation of the operator of Eq.|T|), involving spin 
and isospin degrees of freedom, extremely computation- 
ally intensive, is thus beyond the aim of the present con- 
tribution and will be performed elsewhere. In the present 
work we randomly assigned isospin degrees of freedom for 
A nucleons, so that each generated configuration consists 
of A quads of numbers, containing three spatial coordi- 



Figure 2: The probabilities Pjv, N = 1,2, 3, as a function of the 
impact parameter b (see text), calculated for y/s = 42 GeV (top 
panels) and for y/s = 5500 GeV (LHC energies; lower panels) on the 
ls O nucleus. Squares: analytic results obtained with the approxima- 
tion of Eq. l|12| l. the Glauber approach with zero-range interaction of 
Bertocchi and Treleani, Ref. |15| ; dashes: uncorrelated result; solid 
lines: Gaussian correlations. 



nates and one, randomly assigned, isospin variable; the 
correlated nuclear configurations are downloadable from 
the PSU physics department web server at the address [l3j . 

3. Probability distribution functions for hadron- 
nucleus collisions 

In order to illustrate the present implementation of cor- 
relations, we have focused on the le O, 40 Ca and 20S Pb 
nuclei. 

For a given configuration of nucleons, the probability of 
interaction with the i-th nucleon for an incoming projectile 
with impact parameter b is given by 



P(b,bi) = !-[!- r(6-b l )] i 



(7) 



the corresponding probability of no interaction being 1 — 
P(b, bi). The T function in Eq.([7]) for high energy incident 
nucleons can be parameterized as 
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where cr^^ is the total cross section of NN scattering and 
the t dependence of the cross section da/dt oc exp(Bt), 
neglecting small corrections due to the real part of the 
amplitude. Let us define the probability of interaction 
with N nucleons as a function of the impact parameter as 



Pn (b) = 



N 



A-N 



p(b,b lN ) n [i - p^)] . (9) 

j^i!,...,i N 



0.4 



208 




Figure 3: The probabilities Pjv, N = 1,3,6,9,12, as a function of 
the impact parameter b, calculated for ^/s = 5500 GeV, correspond- 
ing LHC energies, on the 208 Pb nucleus, using (gaussian) correlated 
configurations. 



The inelastic cross section due to collisions with A nucle- 
ons is then given by 



' N 



dbP N {b), 



and the total inelastic cross section is a 
and it can also be calculated as 



NA 



(10) 



'NA 



'NA 



'NA 



'NA ' 



(11) 



0% A an( l a NA being the elastic and quasi-elastic cross sec- 
tion, respectively; this calculation is being performed with 
realistic wave functions (Ref.[3|) and a comparison with 
the present approach will be presented elsewhere. 

We can evaluate the probabilities given by Eq. @ to any 
order in A < A as a function of b = \b\. The results for 
A = 16, 208 are shown in FigsfJl 03 The number of con- 
figurations used in all the calculations is large enough to 
produce negligible statistical errors in the results, typically 
over 100 thousands configurations. All the results have 
been obtained using parameters of the NN amplitude cor- 
responding to = 42 GeV/c (a to t — 41.6 mb, B = 12.6 
GeV~ 2 ), except the lower panel of Fig C3 for 16 0, where 
y/a = 5500 GeV/c (cr tot = 94.8 mb, B = 17.3 GeV- 2 ) is 
chosen, corresponding to LHC energies, to be compared 
with the top panel of the same figure, and in FigJ3] for 
208 Pb where = 5500 GeV/c as well. 

We have checked the independent particle model re- 
sults with the analytic approximation of Ref . . where 
the probability distributions are given by: 



P N (b) 



A\ 



(A-N)\N\ 



(a m T(b)) N [l - a m T{b)] 1 



(12) 



where T(b) = dzp(b,z), which neglects the finite 

radius of NN interaction. The comparison is made for the 
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Figure 4: The dispersion D(b), as a function of the impact parameter 
b, defined in Eq.(Q3 calculated for 16 O, 40 Ca and 208 Pb nuclei, and 
for a 920 GeV incident nucleon. Dashed curves correspond to uncor- 
related configurations; solid lines to the Gaussian correlations. The 
results obtained using eikonal model expression of B&T (Ref. [l5ll ) for 
the zero-range NN interactions are given by dotted lines; for large 
A they coincide with the uncorrelated MC calculations. Symbols: 
direct calculation using Eqs. Ijl3| t . 1140 . from Ref. Q, with C(r) ex- 
tracted from the MC configurations and shown in Fig[T] 



le O nucleus in FigCfl showing a more narrow distribution 
if the approximation of Eq. (fT2"|) is used; the same holds 
true for the other considered nuclei. The quantities (N) 
and (N(N — 1)), can be written as follows[6j: 



(JV> 



drip(ri) 



da" 
~db 



■(6-6i), 



(13) 



(A(A-l)) 



d ri dr 2 p( 2 )(n, r 2 )— -(6 - bx)— -(b - 6a) (14) 
do do 

which can be calculated in the present framework using 



da 1 



db 



(6-6;) = i - (i - r(b-bi)y 



(15) 



Alternatively we can calculate (N), (N(N — 1)) using 
Eqs. (HI), Q3JI as 



N 

(N(N — 1)) = ]T (A 2 - A) P N (b) 



(16) 
(17) 
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Figure 5: The probability density P(x) for the quantity x defined 
in Eq, l|24| l calculated as a function of the impact parameter b, for 



mj = 0.6 GeV 2 and for 208 Pb. We have used the dipole (top) and 



Gaussian (left) gluon density distributions of Ref. [17] and correlated 
configurations. 

As a result we can write the variance of the mean number 
of collisions as a function of the impact parameter b, as 
follows 



D(b) 



(n 2 ) - [(Apr 

(TV) 

E N N 2 P N (b) - [E N NP N (b)} 
En n p nQ>) 



(18) 



In order to check the accuracy of our method, D(b) was 
calculated both using Eqs. (flll]) . (TT4f and Eq. (fT5]) using an- 
alytical one-body densities and C(r) extracted from the 
MC calculation; we then used Eqs. (fTtj|) . (fT7|) . calculated 
with the P/v(6) functions obtained with the MC. The re- 
sults are compared in FigQ] for 16 0, 40 Ca and 208 Pb; a 
small discrepancy is exhibited, between the calculation 
with C(r) (shown with symbols) and the corresponding 
Gaussian correlated result, but the overall agreement is 
satisfactory. The small difference between the results of 
the two methods, especially for lead, is to be ascribed to 
the fact that the quantity C(r) has been extracted from 
the MC configurations using Eq. © , assuming the function 
C( r i: r 2) to be a function of the relative distance only, 
and not to depend on the particular region of the nucleus 
considered, while surface effects should be present on the 
level of few %. We have compared the performance of the 
Glauber type interaction with the approach of Ref. (l6j . 



1.00 
^0.75 
^ 0.50 
0.25 
0.00 



16 °/ 


i 1 i 1 i 1 l 1 




208 Pb y - 




Gaussian 




Dipole 



012345678 

b [fm] 



Figure 6: The ratio R g (b) defined in Eq.(25j for 16 O and 208 Pb, 
calculated using the gaussian (dots) and dipole (solid) gluon density 
distributions of Ref. [Tfl | and correlated configurations as a function 
of the impact parameter b, plotted for m 2 = 0.6 GeV 2 and for those 
values of b for which there is enough statistics for the normalization 
f dxP(x) not to deviate from unity for more than 5%. 



where the probability of interaction in Eq.© was taken 
as a step function, which vanishes for nucleons sitting out- 
side the cylinder centered at the given |b| from the center 
of the nucleus. We took for the cylinder the value which 
gives an area corresponding to the one given by 



dba m (b) 



db 



1 



(i-nb-b,)) 2 (19) 



where T is calculated with the Glauber parameters we have 
used in the present work. The results, as compared with 
our Glauber approach showed to differ substantially from 
both the uncorrelated and correlated MC results. 



4. Lower limit on the parton nuclear shadowing 

In this section we present application of our MC code to 
the study of maximal possible nuclear shadowing in nuclei. 
We use the well known observation that shadowing in the 
scattering off the deuteron cannot reduce the cross section 
to the value smaller than cross section of scattering off 
one nuclcon. Basically it is due to one nucleon screening 
another one but not itself. 

Similarly, it is natural to expect that in any dynam- 
ics of parton interactions the gluon (quark) density at a 
given impact parameter in a particular configuration can- 
not be less than the maximum of individual transverse 
gluon densities, gN(x,p) in the nucleons at given b. Here 
gN(x,p) = gN{x)Fg(p) is the generalized diagonal gluon 
density in the nucleon. 



g A (x, b) co „/ > max i=1 ^A{g N ( x i r t ~ b)) . 



(20) 
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Figure 8: The quantity < x > as a function of A, calculated with 
Figure 7: The same quantity of Fig. [6] but plotted at b = and as a th e densities of Eqs.Q and large values of A. The values of m\ are 
function of m,g, which is a parameter in Eq, l|25|l . given in GeV 2 . 



leading to 

9A(x,b) min > (g A (x,b)conf) , 

In the limit of A — > oo this leads to 
g A (x,b) > g N (p = 0). 



(21) 



(22) 



The onset of the limiting behavior depends on the trans- 
verse shape of the gluon GPD. For the same gN(x) a shape 
more peaked at small p will asymptotically lead to larger 
value of g A (x,b) mm though the approach to the asymp- 
totic value will require larger values of A as chances that 
there is a nucleon with p small enough that Qn(x, p) is close 
to <7at(cc, 0) is smaller in this case. In Ref. [17l| . two param- 
eterizations for the two-gluon form factor were discussed, 
fitted to the J ftp photoproduction data ([3]); they were 
taken in the forms of an exponential and a dipole. The cor- 
responding transverse spatial distributions of gluon GPDs 
are thus 
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where K\ is the modified Bessel function, m 2 — 0.6 GeV 2 
x ~ 10~ 4 and B g = 3.24/m 2 . Using the configurations 



for 



described in the previous sections, we have calculated for 
given impact parameter b the maximum transverse gluon 
density normalized to its peak value, i.e. the maximum 
value of the probability density P(x, b) of the quantity 



F g ( Pi )/F g (0), 



(24) 



as a function of = b — bj, with bj the j-th nucleon trans- 
verse coordinate; here P(x, b) is normalized according to 



J dxP(b, x) = 1, for given b. The results of the calculation 
for 20S Pb are shown in FigEl It can be seen that the two 
densities of Eq. (|23|) produce very similar results. 

We next apply these results to determine minimal value 
of the gluon shadowing at given b defined as the ratio of 
lower limit on g A {x, b) and its value in the impulse approx- 
imation, g A {x,b) — gj^(x)TA(b); to this end, we take the 
ratio 



R g (b) 



<7at(x,p = 0) JdxxP(x,b) 

g N (x)T A (b) 
F g (0) J dxxP(x,b) 



(25) 



Results are presented for the two used models of the gluon 
GPDs of Eq.® in FigEland FigEl For A -> oo the lim- 
its using the two different fits differ by i^ (1) (0)/F 9 (2) (0) = 
2/3.24 = 0.62; however, this limit is reached at extremely 
large A, as shown in FigJS] This is due to a rather 
small radius of the transverse gluon density r* r < 0.5 
fm and low nuclear density, leading to a small probabil- 
ity for more than three nucleons to significantly screen 
each other up to very large A. To investigate the onset of 
asymptotic, we evaluated the A dependence of the quan- 
tity (x) = J dxP{x) at zero impact parameter for the two 
considered gluon GPDs for two values of m 2 correspond- 
ing to x ~ 0.01 and 0.6 GeV 2 corresponding toi~ 10~ 4 . 
We modeled nuclei with very large A by generating ran- 
dom nucleons with constant density po = 0.17 /m -3 in 
a cylinder centered at b = 0, with height 2R A , where 
R A — 1.25A 1 / 3 ; the results, shown in FigEl are indepen- 
dent from the radius of the cylinder, provided it is larger 
than about 0.8 fm. It can be seen that, the GPDs being 
very peaked at p = 0, the increase of (x) is very slow and 
the maximum value is not reached even with A as large as 



6 



10 5 . Onset of asymptotic is somewhat faster for smaller 
m g due to a smoother behavior of GPDs at p ~ 0. 

5. Conclusions 

We developed a MC event generator for nucleon configu- 
rations in nuclei which correctly reproduces single nucleon 
densities and central NN correlations in nuclei. The gener- 
ator can be used in modeling a wide range of processes. In 
particular, it would be interesting to explore how it would 
modify effects of fluctuations of the number of wounded 
nucleons in the heavy ion collisions. 
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